"""
分布模块：药物分布过程建模

该模块实现了药物在体内的分布过程，包括：
- 组织-血浆分配系数计算
- 分布容积计算
- 组织浓度时间曲线
- 生理药代动力学分布模型

作者：BlackCat@CPPO
版本：1.0.0
"""

import numpy as np
from scipy.integrate import solve_ivp
import warnings
import math

class DistributionModel:
    """
    药物分布模型类

    该类模拟药物在体内的分布过程，包括各组织的浓度变化。
    """

    def __init__(self, parameters):
        """
        初始化分布模型

        参数:
            parameters (dict): 模型参数，包含：
                - Vp: 血浆容积 (L)
                - Ve: 间质液容积 (L)
                - Vt: 组织容积 (L)
                - fu: 游离分数
                - BP: 血液-血浆浓度比
                - Kp: 组织-血浆分配系数
                - Q: 组织血流量 (L/h)
        """
        self.params = parameters
        self.tissue_names = [
            'adipose', 'bone', 'brain', 'gut', 'heart', 'kidney',
            'liver', 'lung', 'muscle', 'skin', 'spleen', 'rest'
        ]
        self.validate_parameters()

    def validate_parameters(self):
        """验证模型参数的有效性"""
        required_params = ['Vp', 'Ve', 'fu', 'BP']
        for param in required_params:
            if param not in self.params:
                raise ValueError(f"缺少必需参数：{param}")

    def calculate_distribution(self, individual_params, time_points, absorption_result):
        """
        计算药物分布过程

        参数:
            individual_params (dict): 个体特定参数
            time_points (np.array): 时间点数组
            absorption_result (dict): 吸收阶段结果

        返回:
            dict: 分布过程结果
        """
        # 合并个体参数
        params = {**self.params, **individual_params}

        # 计算组织-血浆分配系数
        tissue_partition = self.calculate_tissue_partition_coefficients(params)

        # 计算分布容积
        distribution_volume = self.calculate_V_d(params, tissue_partition)

        # 计算组织浓度
        tissue_concentrations = self.calculate_tissue_concentrations(
            params, time_points, absorption_result, tissue_partition
        )

        # 计算血浆浓度
        plasma_concentration = self.calculate_plasma_concentration(
            params, time_points, absorption_result, tissue_partition
        )

        # 计算脑部分布
        if isinstance(absorption_result, dict) and 'plasma_concentration' in absorption_result:
            plasma_conc = absorption_result['plasma_concentration']
        else:
            plasma_conc = plasma_concentration

        brain_result = self.calculate_brain_distribution(
            params, time_points, plasma_conc
        )
        tissue_concentrations['brain'] = brain_result['concentration']
        tissue_partition['brain'] = brain_result['partition_coefficient']

        return {
            'tissue_partition': tissue_partition,
            'distribution_volume': distribution_volume,
            'tissue_concentrations': tissue_concentrations,
            'plasma_concentration': plasma_concentration,
            'brain_concentration': brain_result['concentration'],
            'brain_partition': brain_result['partition_coefficient']
        }

    def calculate_tissue_partition_coefficients(self, params):
        """
        计算组织-血浆分配系数

        参数:
            params (dict): 模型参数

        返回:
            dict: 各组织的分配系数
        """
        # 基本参数
        fu = params['fu']
        BP = params['BP']

        # 组织组成参数（体积分数）
        tissue_composition = self._get_tissue_composition()

        # 计算各组织的分配系数
        tissue_partition = {}

        for tissue in self.tissue_names:
            if tissue in tissue_composition:
                comp = tissue_composition[tissue]
                Kp = self._calculate_kp_single_tissue(
                    fu, BP, comp, params
                )
                tissue_partition[tissue] = Kp
            else:
                # 使用默认值
                tissue_partition[tissue] = 1.0

        return tissue_partition

    def _get_tissue_composition(self):
        """获取组织组成数据"""
        return {
            'adipose': {'neutral_lipid': 0.79, 'phospholipid': 0.02, 'water': 0.19},
            'bone': {'neutral_lipid': 0.05, 'phospholipid': 0.01, 'water': 0.44},
            'brain': {'neutral_lipid': 0.12, 'phospholipid': 0.08, 'water': 0.80},
            'gut': {'neutral_lipid': 0.04, 'phospholipid': 0.03, 'water': 0.78},
            'heart': {'neutral_lipid': 0.02, 'phospholipid': 0.01, 'water': 0.77},
            'kidney': {'neutral_lipid': 0.04, 'phospholipid': 0.04, 'water': 0.77},
            'liver': {'neutral_lipid': 0.05, 'phospholipid': 0.05, 'water': 0.71},
            'lung': {'neutral_lipid': 0.01, 'phospholipid': 0.04, 'water': 0.79},
            'muscle': {'neutral_lipid': 0.01, 'phospholipid': 0.01, 'water': 0.76},
            'skin': {'neutral_lipid': 0.02, 'phospholipid': 0.01, 'water': 0.72},
            'spleen': {'neutral_lipid': 0.01, 'phospholipid': 0.05, 'water': 0.78},
            'rest': {'neutral_lipid': 0.03, 'phospholipid': 0.02, 'water': 0.75}
        }

    def _calculate_kp_single_tissue(self, fu, BP, composition, params):
        logP = params.get('logP', 2.0)
        pKa = params.get('pKa', 7.4)
        pH = 7.4  # 生理 pH

        # 计算非离子化分数
        if pKa > pH:
            # 酸性药物：非离子化为主
            non_ionized = 1 / (1 + 10 ** (pH - pKa))
        else:
            # 碱性药物：离子化为主
            non_ionized = 1 / (1 + 10 ** (pKa - pH))

        # 计算分布系数 logD
        logD = logP - math.log10(1 + 10 ** (pH - pKa)) if pKa > pH else logP - math.log10(1 + 10 ** (pKa - pH))

        # 使用 logD 计算各组分贡献
        neutral_lipid_contrib = composition['neutral_lipid'] * non_ionized * (10 ** logD)
        phospholipid_contrib = composition['phospholipid'] * non_ionized * (10 ** logD) * 0.3
        water_contrib = composition['water'] * (1 - non_ionized)

        numerator = neutral_lipid_contrib + phospholipid_contrib + water_contrib
        denominator = 0.04 * (10 ** logD) + 0.022 * (10 ** logD) * 0.3 + 0.94  # 血浆组成

        Kp = (numerator / denominator) * fu / BP

        return max(Kp, 0.1)  # 防止 Kp 过小

    def calculate_V_d(self, params, tissue_partition):
        """
        计算稳态分布容积

        参数:
            params (dict): 模型参数
            tissue_partition (dict): 组织分配系数

        返回:
            float: 稳态分布容积 (L/kg)
        """
        # 基本容积
        Vp = params['Vp']  # 血浆容积
        Ve = params['Ve']  # 间质液容积
        BP = params['BP']

        # 获取组织容积（L/kg体重）
        tissue_volumes = self._get_tissue_volumes()

        # 计算分布容积
        V_d = Vp + Ve

        for tissue, volume in tissue_volumes.items():
            if tissue in tissue_partition:
                V_d += volume * tissue_partition[tissue] / BP

        return V_d

    def _get_tissue_volumes(self):
        """获取组织容积（L/kg体重）"""
        return {
            'adipose': 0.21,
            'bone': 0.05,
            'brain': 0.02,
            'gut': 0.017,
            'heart': 0.005,
            'kidney': 0.004,
            'liver': 0.026,
            'lung': 0.007,
            'muscle': 0.40,
            'skin': 0.037,
            'spleen': 0.002,
            'rest': 0.20
        }

    def calculate_tissue_concentrations(self, params, time_points, absorption_result, tissue_partition):
        """
        计算组织浓度时间曲线

        使用多室PBPK模型计算药物在各组织中的浓度随时间变化

        参数:
            params (dict): 模型参数
            time_points (np.array): 时间点数组
            absorption_result (dict): 吸收结果，包含剂量和时间信息
            tissue_partition (dict): 组织分配系数

        返回:
            dict: 组织浓度时间曲线，包含每个组织的浓度数组
        """
        # 验证输入参数
        if not isinstance(time_points, np.ndarray):
            time_points = np.array(time_points)

        # 获取基本参数
        blood_flows = self._get_blood_flows()
        tissue_volumes = self._get_tissue_volumes()
        fu = params.get('fu', 0.1)
        BP = params.get('BP', 1.0)

        # 获取剂量信息
        if isinstance(absorption_result, dict):
            if 'dose' in absorption_result:
                dose = absorption_result['dose']
            elif 'cumulative' in absorption_result:
                dose = np.max(absorption_result['cumulative'])
            elif 'plasma_concentration' in absorption_result:
                # 如果提供了血浆浓度，使用它来反推剂量
                max_plasma = np.max(absorption_result['plasma_concentration'])
                Vp = params.get('Vp', 3.0)
                dose = max_plasma * Vp * 0.1  # 简化反推
            else:
                dose = 100.0  # 默认剂量
        else:
            # 如果absorption_result是数组，使用它作为输入
            dose = 100.0

        # 初始化结果字典
        tissue_concentrations = {}

        # 计算血浆容积
        plasma_volume = params.get('Vp', 3.0)  # L

        # 对每个组织计算浓度时间曲线
        for tissue in self.tissue_names:
            if tissue in tissue_partition and tissue in blood_flows and tissue in tissue_volumes:
                Kp = tissue_partition[tissue]
                Qt = blood_flows[tissue]  # 组织血流量
                Vt = tissue_volumes[tissue]  # 组织容积

                # 计算组织浓度
                tissue_concentrations[tissue] = self._calculate_single_tissue_concentration(
                    time_points, dose, Kp, Qt, Vt, plasma_volume, fu, BP
                )
            else:
                # 对于缺失参数的组织，使用简化模型
                tissue_concentrations[tissue] = np.zeros_like(time_points)

        return tissue_concentrations

    def _calculate_single_tissue_concentration(self, time_points, dose, Kp, Qt, Vt, Vp, fu, BP):
        """
        计算单个组织的浓度时间曲线

        使用简化但有效的分布模型
        """
        # 计算分配速率
        k_in = (Qt * fu) / (Vt * BP)  # 进入组织的速率
        k_out = Qt / (Vt * Kp)        # 离开组织的速率

        # 计算血浆浓度（简化模型）
        plasma_conc = (dose / Vp) * np.exp(-0.1 * time_points)  # 简化消除

        # 计算组织浓度
        # 使用组织分配平衡模型
        if np.isclose(k_in, k_out):
            tissue_conc = plasma_conc * Kp * (1 - np.exp(-k_in * time_points))
        else:
            tissue_conc = (k_in / (k_out - k_in)) * (
                    plasma_conc * (np.exp(-k_in * time_points) - np.exp(-k_out * time_points))
            )

        # 确保非负值
        tissue_conc = np.maximum(tissue_conc, 0.0)

        return tissue_conc

    def _get_blood_flows(self):
        """获取组织血流量（L/h/kg体重）"""
        return {
            'adipose': 0.026,
            'bone': 0.02,
            'brain': 0.05,
            'gut': 0.15,
            'heart': 0.084,
            'kidney': 0.31,
            'liver': 0.23,
            'lung': 2.4,  # 心输出量
            'muscle': 0.075,
            'skin': 0.026,
            'spleen': 0.056,
            'rest': 0.03
        }

    def calculate_plasma_concentration(self, params, time_points, absorption_result, tissue_partition):
        """
        计算血浆浓度

        参数:
            params (dict): 模型参数
            time_points (np.array): 时间点数组
            absorption_result (dict): 吸收结果
            tissue_partition (dict): 组织分配系数

        返回:
            np.array: 血浆浓度时间曲线
        """
        # 计算分布容积
        V_d = self.calculate_V_d(params, tissue_partition)

        # 获取吸收参数
        ka = params.get('ka', 0.5)  # 吸收速率常数

        # 从吸收结果中获取数据
        if isinstance(absorption_result, dict) and 'cumulative' in absorption_result:
            cumulative_absorption = absorption_result['cumulative']
            dose = params.get('dose', 100.0)
        else:
            # 如果absorption_result是数组，使用它作为累积吸收量
            cumulative_absorption = np.array(absorption_result)
            dose = np.max(cumulative_absorption) if cumulative_absorption.size > 0 else 100.0

        # 使用一室模型计算血浆浓度
        # 计算消除速率常数（基于分布容积和清除率）
        CL_total = params.get('CL_hepatic', 0.5) + params.get('CL_renal', 0.3) + params.get('CL_biliary', 0.1)
        ke = CL_total / V_d  # 消除速率常数

        # 确保ka和ke不相等
        if abs(ka - ke) < 1e-6:
            ke = ke * 1.01  # 稍微调整ke以避免除零

        # 使用一室模型公式计算血浆浓度
        # C(t) = (F * D * ka) / (V * (ka - ke)) * (e^(-ke*t) - e^(-ka*t))
        # 其中F是生物利用度
        if isinstance(absorption_result, dict) and 'bioavailability' in absorption_result:
            F = absorption_result['bioavailability']
        else:
            F = params.get('fa', 0.8) * params.get('fg', 0.9) * params.get('fh', 0.7)

        # 处理吸收与消除速率接近的情况
        if np.isclose(ka, ke, atol=1e-3):
            # 简化为单指数衰减（静脉注射情况）
            plasma_concentration = (dose / V_d) * np.exp(-ke * time_points)
        else:
            # 计算血浆浓度
            plasma_concentration = (F * dose * ka) / (V_d * (ka - ke)) * (
                np.exp(-ke * time_points) - np.exp(-ka * time_points)
            )

        # 确保浓度非负
        plasma_concentration = np.maximum(plasma_concentration, 0.0)

        return plasma_concentration

    def calculate_tissue_distribution_time(self, tissue):
        """
        计算组织分布时间

        参数:
            tissue (str): 组织名称

        返回:
            float: 分布半衰期 (小时)
        """
        blood_flows = self._get_blood_flows()
        tissue_volumes = self._get_tissue_volumes()

        if tissue not in blood_flows or tissue not in tissue_volumes:
            return np.nan

        # 分布半衰期 = 0.693 * 组织容积 / 血流量
        t_half = 0.693 * tissue_volumes[tissue] / blood_flows[tissue]

        return t_half

    def calculate_brain_distribution(self, params, time_points, plasma_concentration):
        """
        计算脑部分布（特殊考虑血脑屏障）

        参数:
            params (dict): 模型参数
            time_points (np.array): 时间点数组
            plasma_concentration (np.array): 血浆浓度

        返回:
            dict: 脑部浓度结果
        """
        # 血脑屏障参数
        bbb_permeability = params.get('bbb_permeability', 0.01)  # cm/s
        brain_flow = 0.05  # L/h/kg
        brain_volume = 0.02  # L/kg

        # 计算脑部分配系数（考虑血脑屏障）
        Kp_brain = self._calculate_brain_partition(params)

        # 计算脑部浓度（简化模型）
        brain_concentration = plasma_concentration * Kp_brain * (
            1 - np.exp(-bbb_permeability * brain_flow / brain_volume * time_points)
        )

        return {
            'concentration': brain_concentration,
            'partition_coefficient': Kp_brain,
            'equilibrium_time': 3.0 / (bbb_permeability * brain_flow / brain_volume)
        }

    def _calculate_brain_partition(self, params):
        """计算脑部分配系数（考虑 pKa 和 pH）"""
        logP = params.get('logP', 2.0)
        pKa = params.get('pKa', 7.4)
        fu = params.get('fu', 0.1)  # 游离分数
        pH = 7.4  # 生理 pH

        # 计算非离子化比例（假设为弱酸性药物）
        if pKa > pH:
            non_ionized_fraction = 10 ** (pKa - pH)
        else:
            non_ionized_fraction = 1 / (10 ** (pH - pKa))

        # 计算 logD（分布系数）
        logD = logP - math.log10(1 + 10 ** (pH - pKa)) if pKa > pH else logP - math.log10(1 + 10 ** (pKa - pH))

        # 计算血液-脑分配系数（P_blood_brain）
        P_blood_brain = 0.001 * (logD ** 3) - 0.03 * (logD ** 2) + 0.15 * logD + 0.02

        # 计算脑部分配系数
        Kp_brain = (fu * non_ionized_fraction) * P_blood_brain

        # 添加阈值限制（可选）
        Kp_brain = max(Kp_brain, 0.01)  # 最小值设为 0.01
        
        return Kp_brain
    
    def calculate_fetal_distribution(self, params, maternal_concentration):
        """
        计算胎儿分布（孕期用药）

        参数:
            params (dict): 模型参数
            maternal_concentration (np.array): 母体血浆浓度 (mg/L)

        返回:
            dict: 胎儿浓度结果
        """
        # 胎盘屏障参数
        placental_permeability = params.get('placental_permeability', 0.1)  # cm/s
        fetal_blood_flow = params.get('fetal_blood_flow', 0.15)  # L/h/kg
        fetal_volume = params.get('fetal_volume', 0.1)  # L/kg
        fetal_elimination_rate = params.get('fetal_elimination_rate', 0.01)  # 1/h

        # 单位转换：cm/s → L/(h·kg)
        placental_permeability = placental_permeability * 3600 * 1e-6  # L/(h·kg)

        # 计算胎盘转移率 (L/h/kg)
        transfer_rate = placental_permeability * fetal_blood_flow

        # 计算胎儿浓度（动态平衡模型）
        # dC/dt = (transfer_rate * maternal_concentration) - (fetal_elimination_rate * C)
        # 解析解：C(t) = (transfer_rate / (fetal_elimination_rate)) * maternal_concentration
        # 假设稳态或快速平衡（实际需用 ODE 求解器更精确）
        fetal_concentration = (transfer_rate * maternal_concentration) / (fetal_elimination_rate + transfer_rate)

        return {
            'concentration': fetal_concentration,
            'fetal_volume': fetal_volume,
            'fetal_blood_flow': fetal_blood_flow,
            'placental_transfer': placental_permeability,
            'fm_ratio': fetal_elimination_rate
        }

    def calculate_multi_compartment_plasma_concentration(
            self,
            params: dict,  # 包含生理参数（如 Vp, Q_liver, CL_hepatic）和可选的 body_weight
            time_points: np.array,
            absorption_result: dict
    ) -> dict:
        """
        使用完整多室PBPK模型计算血浆和组织浓度
        
        实现包含主要器官和组织的完整生理药代动力学模型：
        - 血浆（中央室）
        - 肝脏（代谢器官）
        - 肾脏（排泄器官）
        - 脑、心脏、肺、肌肉、脂肪、骨骼等
        
        参数:
        params (dict): 必须包含：
            - Vp: 血浆容积 (L/kg)
            - Q_liver: 肝血流 (L/h/kg)
            - CL_hepatic: 肝清除率 (L/h)
            - ...其他器官参数...
            - body_weight: 个体体重 (kg, 可选，默认70)
        time_points (np.array): 时间点数组
        absorption_result (dict): 吸收阶段结果（剂量、生物利用度等）
            
        返回:
            dict: 完整多室模型结果，包含血浆浓度和各组织浓度
        """
        
        # 获取基本生理参数
        body_weight = params.get('body_weight', 70.0)  # kg
        
        # 定义组织容积（L/kg体重）- 调整为与最小模型一致
        # 使用实际生理容积，但确保血浆容积与最小模型一致
        plasma_vol = params.get('Vp', 3.0) / body_weight # 转换为L/kg
        tissue_volumes = {
            'plasma': plasma_vol,
            'liver': params.get('V_liver', 0.026),
            'kidney': params.get('V_kidney', 0.004),
            'brain': params.get('V_brain', 0.02),
            'heart': params.get('V_heart', 0.005),
            'lung': params.get('V_lung', 0.007),
            'muscle': params.get('V_muscle', 0.40),
            'adipose': params.get('V_adipose', 0.21),
            'bone': params.get('V_bone', 0.05),
            'gut': params.get('V_gut', 0.017),
            'spleen': params.get('V_spleen', 0.0026),
            'skin': params.get('V_skin', 0.037),
            'rest': params.get('V_rest', 0.10)
        }
        
        # 定义血流量（L/h/kg）
        cardiac_output = params.get('cardiac_output', 5.0)  # L/h/kg
        blood_flows = {
            'liver': params.get('Q_liver', 0.25 * cardiac_output),  # 肝脏血流量
            'kidney': params.get('Q_kidney', 0.19 * cardiac_output),  # 肾脏血流量
            'brain': params.get('Q_brain', 0.12 * cardiac_output),   # 脑血流量
            'heart': params.get('Q_heart', 0.04 * cardiac_output),   # 心脏血流量
            'lung': params.get('Q_lung', 1.0 * cardiac_output),      # 肺血流量（100%）
            'muscle': params.get('Q_muscle', 0.17 * cardiac_output), # 肌肉血流量
            'adipose': params.get('Q_adipose', 0.05 * cardiac_output), # 脂肪血流量
            'bone': params.get('Q_bone', 0.05 * cardiac_output),     # 骨骼血流量
            'gut': params.get('Q_gut', 0.17 * cardiac_output),       # 肠道血流量
            'spleen': params.get('Q_spleen', 0.03 * cardiac_output), # 脾脏血流量
            'skin': params.get('Q_skin', 0.05 * cardiac_output),     # 皮肤血流量
            'rest': params.get('Q_rest', 0.12 * cardiac_output)      # 其余组织血流量
        }
        
        # 计算实际容积（L）- 使用直接容积而非按体重缩放
        actual_volumes = {tissue: vol * body_weight for tissue, vol in tissue_volumes.items()}  # 固定70kg成人体重
        actual_flows = {tissue: flow * body_weight for tissue, flow in blood_flows.items()}     # 固定70kg成人体重
        
        # 获取分配系数
        tissue_partition = self.calculate_tissue_partition_coefficients(params)
        
        # 使用与最小模型一致的清除率（已按70kg成人体重标准化）
        CL_hepatic = max(params.get('CL_hepatic', 2.0), 0.1)  # L/h
        CL_renal = max(params.get('CL_renal', 1.5), 0.1)      # L/h
        CL_biliary = max(params.get('CL_biliary', 0.5), 0.05)  # L/h
        
        # 获取剂量和吸收参数
        dose = params.get('dose', 100.0)
        ka = params.get('ka', 0.5)
        
        # 获取生物利用度
        if isinstance(absorption_result, dict) and 'bioavailability' in absorption_result:
            F = absorption_result['bioavailability']
        else:
            F = params.get('fa', 0.8) * params.get('fg', 0.9) * params.get('fh', 0.7)
        
        # 定义组织列表
        tissues = ['plasma', 'liver', 'kidney', 'brain', 'heart', 'lung', 
                  'muscle', 'adipose', 'bone', 'gut', 'spleen', 'skin', 'rest']
        
        # 定义状态变量：每个组织的药物量
        def pbpk_model(y, t):
            """
            完整多室PBPK模型的微分方程组
            y: [A_plasma, A_liver, A_kidney, A_brain, A_heart, A_lung, A_muscle, 
                A_adipose, A_bone, A_gut, A_spleen, A_skin, A_rest]
            其中A是药物量(mg)，不是浓度(mg/L)
            """
            # 解包状态变量（药物量），添加边界保护
            amounts = {}
            for i, tissue in enumerate(tissues):
                amounts[tissue] = max(y[i], 1e-12)  # 防止数值问题
            
            # 计算各组织的浓度，确保非负
            concentrations = {}
            for tissue in tissues:
                concentrations[tissue] = max(amounts[tissue] / actual_volumes[tissue], 1e-12)
            
            # 计算吸收输入（随时间变化），使用更安全的计算
            absorption_t = max(F * dose * ka * np.exp(-ka * t), 1e-12) if t >= 0 else 0
            
            # 计算各组织的清除率（基于组织浓度，单位：mg/h）
            # 注意：这里应该使用血浆浓度作为清除的基础，因为清除主要发生在血浆
            total_clearance = CL_hepatic + CL_renal + CL_biliary
            clearance_plasma = total_clearance * concentrations['plasma']  # L/h * mg/L = mg/h
            
            # 器官特异性清除（简化模型）
            clearance_liver = CL_hepatic * concentrations['liver']  # L/h * mg/L = mg/h
            clearance_kidney = CL_renal * concentrations['kidney']  # L/h * mg/L = mg/h
            clearance_biliary = CL_biliary * concentrations['liver']  # L/h * mg/L = mg/h
            
            # 初始化变化率数组（药物量变化率，mg/h）
            dA_dt = np.zeros(len(tissues))
            
            # 血浆室（中央室）变化率
            dA_dt[0] = absorption_t - clearance_plasma  # 吸收输入 - 系统清除
            
            # 计算组织交换对血浆的影响
            total_exchange_to_plasma = 0
            
            for i, tissue in enumerate(tissues[1:], 1):
                partition_coeff = max(tissue_partition.get(tissue, 1.0), 0.01)
                flow = max(actual_flows.get(tissue, 0), 1e-12)
                
                # 计算浓度梯度
                plasma_conc = concentrations['plasma']
                tissue_conc = concentrations[tissue]
                
                # 血浆→组织交换（基于浓度梯度），确保数值稳定
                concentration_diff = plasma_conc - (tissue_conc / partition_coeff)
                exchange_flux = flow * concentration_diff
                
                # 组织药物量变化
                dA_dt[i] = exchange_flux
                
                # 特定组织的清除处理
                if tissue == 'liver':
                    # 肝脏代谢清除和胆汁排泄清除
                    dA_dt[i] -= (clearance_liver + clearance_biliary)
                elif tissue == 'kidney':
                    # 肾脏排泄清除
                    dA_dt[i] -= clearance_kidney
                
                # 累积对血浆的交换影响
                total_exchange_to_plasma -= exchange_flux
            
            # 添加组织交换对血浆的影响
            dA_dt[0] += total_exchange_to_plasma
            
            # 确保变化率不会过大
            max_rate = 1000.0
            dA_dt = np.clip(dA_dt, -max_rate, max_rate)
            
            return dA_dt
        
        # 初始条件：所有室均为0，药物通过吸收逐步进入系统
        y0 = [0.0] * len(tissues)  # 所有组织初始药物量为0
        
        # 求解微分方程
        try:
            solution = solve_ivp(
                lambda t, y: pbpk_model(y, t), 
                [time_points[0], time_points[-1]], 
                y0, 
                t_eval=time_points,
                rtol=1e-6, 
                atol=1e-6,
                method='RK45'  # 使用Runge-Kutta方法求解
            )
            if solution.success:
                solution = solution.y.T
            else:
                warnings.warn(f"ODE求解失败: {solution.message}")
                return None
        except RuntimeError as e:
            warnings.warn(f"ODE求解失败: {str(e)}")
            return None
        
        # 提取结果（转换为浓度）
        results = {
            'plasma_concentration': solution[:, 0] / actual_volumes['plasma'],
            'tissue_concentrations': {
                tissue: solution[:, i + 1] / actual_volumes[tissue]
                for i, tissue in enumerate(tissues[1:])
            }
        }
        
        # 确保非负值
        results['plasma_concentration'] = np.maximum(results['plasma_concentration'], 0.0)
        for tissue in results['tissue_concentrations']:
            results['tissue_concentrations'][tissue] = np.maximum(
                results['tissue_concentrations'][tissue], 0.0
            )
        
        return results
    
    def _fallback_multi_compartment(self, params, time_points):
        """多室模型求解失败时的备用方案"""
        # 使用二室模型作为备用
        V1 = params.get('Vp', 3.0)  # 中央室容积
        V2 = params.get('Ve', 12.0)  # 周边室容积
        Q = params.get('Q_intercompartmental', 15.0)  # 室间清除率
        CL = params.get('CL_hepatic', 0.5) + params.get('CL_renal', 0.3)  # 总清除率
        
        # 计算速率常数
        k10 = CL / V1
        k12 = Q / V1
        k21 = Q / V2
        
        # 计算混合参数
        lambda1 = 0.5 * ((k10 + k12 + k21) + np.sqrt((k10 + k12 + k21)**2 - 4*k10*k21))
        lambda2 = 0.5 * ((k10 + k12 + k21) - np.sqrt((k10 + k12 + k21)**2 - 4*k10*k21))
        
        # 计算系数
        A = (k10 - lambda2) / (lambda1 - lambda2)
        B = (lambda1 - k10) / (lambda1 - lambda2)
        
        # 获取剂量和吸收参数
        dose = params.get('dose', 100.0)
        ka = params.get('ka', 0.5)
        F = params.get('fa', 0.8) * params.get('fg', 0.9) * params.get('fh', 0.7)
        
        # 计算血浆浓度（二室模型）
        plasma_conc = np.zeros_like(time_points)
        
        for i, t in enumerate(time_points):
            if ka == lambda1 or ka == lambda2:
                # 避免除零
                continue
            
            term1 = A * (np.exp(-lambda1 * t) - np.exp(-ka * t)) / (ka - lambda1)
            term2 = B * (np.exp(-lambda2 * t) - np.exp(-ka * t)) / (ka - lambda2)
            
            plasma_conc[i] = (F * dose * ka / V1) * (term1 + term2)
        
        # 确保非负
        plasma_conc = np.maximum(plasma_conc, 0.0)
        
        # 简化的组织浓度
        tissue_concentrations = {}
        tissue_partition = self.calculate_tissue_partition_coefficients(params)
        
        for tissue in ['liver', 'kidney', 'brain', 'heart', 'muscle', 'adipose']:
            if tissue in tissue_partition:
                # 使用平衡假设计算组织浓度
                tissue_concentrations[tissue] = plasma_conc * tissue_partition[tissue] * (
                    1 - np.exp(-0.1 * time_points)
                )
        
        return {
            'plasma_concentration': plasma_conc,
            'tissue_concentrations': tissue_concentrations
        }